library(dplyr)
library(tidyr)
library(car)
library(readxl)
library(mice)
library(VIM)
library(grid)
library(ggplot2)
library(cowplot)
library(missMDA)
library(FactoMineR)
library(TidyDensity)
library(MASS)
library(leaps)
library(lmtest)
library(Metrics)

library(pls)
# library(caret)

Attaching package: ‘dplyr’


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


Loading required package: carData


Attaching package: ‘car’


The following object is masked from ‘package:dplyr’:

    recode


Attaching package: ‘mice’


The following object is masked from ‘package:stats’:

    filter


The following objects are masked from ‘package:base’:

    cbind, rbind


Loading required package: colorspace

Loading required package: grid

VIM is ready to use.


Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues


Attaching package: ‘VIM’


The following object is masked from ‘package:datasets’:

    sleep


== Welcome to TidyDensity ===========================================================================
If you find this package useful, please leave a star: 
   https://github.com/spsanderson/TidyDensity'

If you encounter a bug or want to request an enhancement please file an issue at:
   https://github.com/spsanderson/TidyDensity/issues

Thank you for using TidyDensity!


Attaching package: ‘MASS’


The following object is masked from ‘package:dplyr’:

    select


Loading required package: zoo


Attaching package: ‘zoo’


The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric


Attaching package: ‘pls’


The following object is masked from ‘package:stats’:

    loadings

remove_high_vif <- function(regr_variable, data, threshold = 5) {      
    stop <- FALSE
    while(stop == FALSE) {
        # Tính VIF cho các biến độc lập
        model <- lm(as.formula(paste(regr_variable, "~ .")), data = data)
        vif_values <- vif(model)
        # Tìm các biến có VIF cao hơn ngưỡng
        high_vif <- which(vif_values > threshold)
        
        if (length(high_vif) > 0) {
            # Loại bỏ các biến có VIF cao
            remove_index <- which.max(high_vif)
            rm_names <- names(data[, !names(data) %in% regr_variable])
            rm_name <- rm_names[remove_index]
            data <- data[, !names(data) %in% rm_name]
        } else {
            stop <- TRUE
        }
    }
    return(data)
}

response_variable = 1 

# Hàm tiền xử lý dữ liệu với box-cox
bc_transform <- function(df) {
    col.names <- names(df)

    transformed_df <- df
    
    for (name in names(df))
    {
        col.name <- name

        print(col.name)

        # Rút trích biến phản hồi
        response_variable <- df[[col.name]]
        
        if (!is.numeric(response_variable)) {
            print(col.name)
            stop("The column to be transformed must be numeric.")
        }
        
        if (any(response_variable <= 0)) {
            # Shift the values to be positive
            shift_value <- abs(min(response_variable)) + 1
            response_variable <- response_variable + shift_value
        }
        
        # Áp dụng box-cox transform để tìm lambda tối ưu
        boxcox_result <- boxcox(lm(response_variable ~ 1), plotit = FALSE)
        optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
        print(paste("Optimal lambda:", optimal_lambda))
        
        # Sử dụng lambda tối ưu để biến đổi dữ liệu
        if (optimal_lambda == 0) {
            transformed_response <- log(response_variable)
        } else {
            transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
        }
        
        # Gán biến đã được biến đổi
        transformed_df[[col.name]] <- transformed_response
    }
    
    return(transformed_df)
}

indicator <- function(model, y_pred, y_true) {
     adj.r.sq <- summary(model)$adj.r.squared
     mse <- MSE(y_pred, y_true)
     rmse <- RMSE(y_pred, y_true)
     mae <- MAE(y_pred, y_true)
     print(paste0("Adjusted R-squared: ", round(adj.r.sq, 4)))
     print(paste0("MSE: ", round(mse, 4)))
     print(paste0("RMSE: ", round(rmse, 4)))
     print(paste0("MAE: ", round(mae, 4)))
}


metrics <- function(y_pred, y_true){
     mse <- MSE(y_pred, y_true)
     rmse <- RMSE(y_pred, y_true)
     mae <- MAE(y_pred, y_true)
     print(paste0("MSE: ", round(mse, 6)))
     print(paste0("RMSE: ", round(rmse, 6)))
     print(paste0("MAE: ", round(mae, 6)))
     corPredAct <- cor(y_pred, y_true)
     print(paste0("Correlation: ", round(corPredAct, 6)))
     print(paste0("R^2 between y_pred & y_true: ", round(corPredAct^2, 6)))
}

CheckNormal <- function(model) {
    if (length(model$residuals) < 4997)
    {
        hist(model$residuals, breaks = 30)
        shaptest <- shapiro.test(model$residuals)
        print(shaptest)
        if (shaptest$p.value <= 0.05) {
            print("H0 rejected: the residuals are NOT distributed normally")
        } else {
            print("H0 failed to reject: the residuals ARE distributed normally")
        }
    } else {
        hist(model$residuals, breaks = 30)
        shaptest <- shapiro.test(model$residuals[3:5000])
        print(shaptest)
        if (shaptest$p.value <= 0.05) {
            print("H0 rejected: the residuals are NOT distributed normally")
        } else {
            print("H0 failed to reject: the residuals ARE distributed normally")
        }
    }
}

CheckHomos <- function(model){
     plot(model$fitted.values, model$residuals)
     abline(h = 0, col = "red")
     BP <- bptest(model)
     print(BP)
     if (BP$p.value <= 0.05) {
          print("H0 rejected: Error variance spreads INCONSTANTLY/generating patterns (Heteroscedasticity)")
     } else {
          print("H0 failed to reject: Error variance spreads CONSTANTLY (Homoscedasticity)")
     }
}

# Đọc dữ liệu từ tập tin
raw_data = read_excel("../../data/part2/air+quality/AirQualityUCI.xlsx", sheet = 1)
str(raw_data)

tibble [9,357 × 15] (S3: tbl_df/tbl/data.frame)
 $ Date         : POSIXct[1:9357], format: "2004-03-10" "2004-03-10" ...
 $ Time         : POSIXct[1:9357], format: "1899-12-31 18:00:00" "1899-12-31 19:00:00" ...
 $ CO(GT)       : num [1:9357] 2.6 2 2.2 2.2 1.6 1.2 1.2 1 0.9 0.6 ...
 $ PT08.S1(CO)  : num [1:9357] 1360 1292 1402 1376 1272 ...
 $ NMHC(GT)     : num [1:9357] 150 112 88 80 51 38 31 31 24 19 ...
 $ C6H6(GT)     : num [1:9357] 11.88 9.4 9 9.23 6.52 ...
 $ PT08.S2(NMHC): num [1:9357] 1046 955 939 948 836 ...
 $ NOx(GT)      : num [1:9357] 166 103 131 172 131 89 62 62 45 -200 ...
 $ PT08.S3(NOx) : num [1:9357] 1056 1174 1140 1092 1205 ...
 $ NO2(GT)      : num [1:9357] 113 92 114 122 116 96 77 76 60 -200 ...
 $ PT08.S4(NO2) : num [1:9357] 1692 1559 1554 1584 1490 ...
 $ PT08.S5(O3)  : num [1:9357] 1268 972 1074 1203 1110 ...
 $ T            : num [1:9357] 13.6 13.3 11.9 11 11.2 ...
 $ RH           : num [1:9357] 48.9 47.7 54 60 59.6 ...
 $ AH           : num [1:9357] 0.758 0.725 0.75 0.787 0.789 ...

# Kích thước của dữ liệu
# Ta thấy dữ liệu có 231 dòng và 14 cột
dim(raw_data)

head(raw_data)

# Kiểm tra dữ liệu trùng lặp
duplicates <- raw_data[duplicated(raw_data), ]
duplicate_counts <- table(raw_data[duplicated(raw_data), ])
duplicates # Không có dữ liệu trùng lặp

# Đầu tiên, ta xem xét một số quan trắc
str(raw_data)

tibble [9,357 × 15] (S3: tbl_df/tbl/data.frame)
 $ Date         : POSIXct[1:9357], format: "2004-03-10" "2004-03-10" ...
 $ Time         : POSIXct[1:9357], format: "1899-12-31 18:00:00" "1899-12-31 19:00:00" ...
 $ CO(GT)       : num [1:9357] 2.6 2 2.2 2.2 1.6 1.2 1.2 1 0.9 0.6 ...
 $ PT08.S1(CO)  : num [1:9357] 1360 1292 1402 1376 1272 ...
 $ NMHC(GT)     : num [1:9357] 150 112 88 80 51 38 31 31 24 19 ...
 $ C6H6(GT)     : num [1:9357] 11.88 9.4 9 9.23 6.52 ...
 $ PT08.S2(NMHC): num [1:9357] 1046 955 939 948 836 ...
 $ NOx(GT)      : num [1:9357] 166 103 131 172 131 89 62 62 45 -200 ...
 $ PT08.S3(NOx) : num [1:9357] 1056 1174 1140 1092 1205 ...
 $ NO2(GT)      : num [1:9357] 113 92 114 122 116 96 77 76 60 -200 ...
 $ PT08.S4(NO2) : num [1:9357] 1692 1559 1554 1584 1490 ...
 $ PT08.S5(O3)  : num [1:9357] 1268 972 1074 1203 1110 ...
 $ T            : num [1:9357] 13.6 13.3 11.9 11 11.2 ...
 $ RH           : num [1:9357] 48.9 47.7 54 60 59.6 ...
 $ AH           : num [1:9357] 0.758 0.725 0.75 0.787 0.789 ...

# 1. Time
is.factor(raw_data$Time) #False 

# 2. PT08.S1(CO)
is.factor(raw_data$`PT08.S1(CO)`) # False

# 3. PT08.S2(NMHC)
is.factor(raw_data$`PT08.S2(NMHC)`) # False

# 4. PT08.S3(NOx)
is.factor(raw_data$`PT08.S3(NOx)`) # False

# 5. PT08.S4(NO2)
is.factor(raw_data$`PT08.S4(NO2)`) # False

# 6. PT08.S5(O3)
is.factor(raw_data$`PT08.S5(O3)`) # False

# Convert character columns to numeric, except for Date and Time
cols_to_convert <- setdiff(names(raw_data), c("Date", "Time"))

# Use mutate to convert columns to numeric
air_quality <- raw_data %>% mutate(across(all_of(cols_to_convert), ~ as.numeric(gsub(",", ".", .))))
head(air_quality)

# # Modifications
# raw_data$Time <- as.factor(raw_data$Time)
# raw_data$`PT08.S1(CO)` <- as.factor(raw_data$`PT08.S1(CO)`)
# raw_data$`PT08.S2(NMHC)` <- as.factor(raw_data$`PT08.S2(NMHC)`)
# raw_data$`PT08.S3(NOx)` <- as.factor(raw_data$`PT08.S3(NOx)`)
# raw_data$`PT08.S4(NO2)` <- as.factor(raw_data$`PT08.S4(NO2)`)
# raw_data$`PT08.S5(O3)` <- as.factor(raw_data$`PT08.S5(O3)`)

# Hàm tính toán tỷ lệ missing values
missing_ratio <- function(s) {
  round(mean(is.na(s)) * 100, 1)
}

# Hàm tính toán trung vị (median)
median_custom <- function(df) {
  round(quantile(df, 0.5, na.rm = TRUE), 1)
}

# Hàm tính toán phân vị 25% (Q1)
lower_quartile <- function(df) {
  round(quantile(df, 0.25, na.rm = TRUE), 1)
}

# Hàm tính toán phân vị 75% (Q3)
upper_quartile <- function(df) {
  round(quantile(df, 0.75, na.rm = TRUE), 1)
}

# Lựa chọn các kiểu số
num_col_info_df <- as.data.frame(air_quality) %>% select_if(is.numeric)

# Tổng hợp thông tin thống kê mô tả
num_col_info_df <- as.data.frame(air_quality) %>%
  select_if(is.numeric) %>%
  summarise(
    across(everything(), list(
      missing_ratio = ~ missing_ratio(.),
      min = ~ min(., na.rm = TRUE),
      lower_quartile = ~ lower_quartile(.),
      median = ~ median_custom(.),
      upper_quartile = ~ upper_quartile(.),
      max = ~ max(., na.rm = TRUE)
    ))
  )

num_col_info_df <- num_col_info_df %>%
  pivot_longer(
    cols = everything(),
    names_to = c("variable", ".value"),
    names_sep = "_"
  )

print(num_col_info_df)

Warning message:
“Expected 2 pieces. Additional pieces discarded in 39 rows [1, 3, 5, 7, 9, 11,
13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, ...].”

# A tibble: 13 × 7
   variable      missing   min  lower median  upper     max
   <chr>           <dbl> <dbl>  <dbl>  <dbl>  <dbl>   <dbl>
 1 CO(GT)              0  -200    0.6    1.5    2.6   11.9 
 2 PT08.S1(CO)         0  -200  921   1052.  1221.  2040.  
 3 NMHC(GT)            0  -200 -200   -200   -200   1189   
 4 C6H6(GT)            0  -200    4      7.9   13.6   63.7 
 5 PT08.S2(NMHC)       0  -200  711    894.  1105.  2214   
 6 NOx(GT)             0  -200   50    141    284.  1479   
 7 PT08.S3(NOx)        0  -200  637    794.   960.  2683.  
 8 NO2(GT)             0  -200   53     96    133    340.  
 9 PT08.S4(NO2)        0  -200 1185.  1446.  1662   2775   
10 PT08.S5(O3)         0  -200  700.   942   1255.  2523.  
11 T                   0  -200   10.9   17.2   24.1   44.6 
12 RH                  0  -200   34.1   48.6   61.9   88.7 
13 AH                  0  -200    0.7    1      1.3    2.23

print(summary(as.data.frame(air_quality) %>% select_if(is.numeric)))

     CO(GT)         PT08.S1(CO)      NMHC(GT)         C6H6(GT)       
 Min.   :-200.00   Min.   :-200   Min.   :-200.0   Min.   :-200.000  
 1st Qu.:   0.60   1st Qu.: 921   1st Qu.:-200.0   1st Qu.:   4.005  
 Median :   1.50   Median :1052   Median :-200.0   Median :   7.887  
 Mean   : -34.21   Mean   :1049   Mean   :-159.1   Mean   :   1.866  
 3rd Qu.:   2.60   3rd Qu.:1221   3rd Qu.:-200.0   3rd Qu.:  13.636  
 Max.   :  11.90   Max.   :2040   Max.   :1189.0   Max.   :  63.741  
 PT08.S2(NMHC)       NOx(GT)        PT08.S3(NOx)       NO2(GT)       
 Min.   :-200.0   Min.   :-200.0   Min.   :-200.0   Min.   :-200.00  
 1st Qu.: 711.0   1st Qu.:  50.0   1st Qu.: 637.0   1st Qu.:  53.00  
 Median : 894.5   Median : 141.0   Median : 794.2   Median :  96.00  
 Mean   : 894.5   Mean   : 168.6   Mean   : 794.9   Mean   :  58.14  
 3rd Qu.:1104.8   3rd Qu.: 284.2   3rd Qu.: 960.2   3rd Qu.: 133.00  
 Max.   :2214.0   Max.   :1479.0   Max.   :2682.8   Max.   : 339.70  
  PT08.S4(NO2)   PT08.S5(O3)           T                  RH         
 Min.   :-200   Min.   :-200.0   Min.   :-200.000   Min.   :-200.00  
 1st Qu.:1185   1st Qu.: 699.8   1st Qu.:  10.950   1st Qu.:  34.05  
 Median :1446   Median : 942.0   Median :  17.200   Median :  48.55  
 Mean   :1391   Mean   : 975.0   Mean   :   9.777   Mean   :  39.48  
 3rd Qu.:1662   3rd Qu.:1255.2   3rd Qu.:  24.075   3rd Qu.:  61.88  
 Max.   :2775   Max.   :2522.8   Max.   :  44.600   Max.   :  88.72  
       AH           
 Min.   :-200.0000  
 1st Qu.:   0.6923  
 Median :   0.9768  
 Mean   :  -6.8376  
 3rd Qu.:   1.2962  
 Max.   :   2.2310

# # Gán các giá trị bé hơn 0, bằng 0
replace_negatives <- function(x) {
  x[x == -200] <- NA
  return(x)
}

air_quality <- as.data.frame(lapply(air_quality, replace_negatives))

tmp <- air_quality[, -c(1, 2)]

# using colMeans()
mean_val <- colMeans(tmp, na.rm = TRUE)
 
# replacing NA with mean value of each column
for(i in colnames(tmp))
  tmp[,i][is.na(tmp[,i])] <- mean_val[i]

tmp$Date <- air_quality$Date
tmp$Time <- air_quality$Time
air_quality <- tmp

print(summary(as.data.frame(air_quality) %>% select_if(is.numeric)))

     CO.GT.        PT08.S1.CO.        NMHC.GT.         C6H6.GT.     
 Min.   : 0.100   Min.   : 647.2   Min.   :   7.0   Min.   : 0.149  
 1st Qu.: 1.200   1st Qu.: 941.2   1st Qu.: 218.8   1st Qu.: 4.591  
 Median : 2.153   Median :1074.5   Median : 218.8   Median : 8.593  
 Mean   : 2.153   Mean   :1099.7   Mean   : 218.8   Mean   :10.083  
 3rd Qu.: 2.600   3rd Qu.:1221.2   3rd Qu.: 218.8   3rd Qu.:13.636  
 Max.   :11.900   Max.   :2039.8   Max.   :1189.0   Max.   :63.741  
 PT08.S2.NMHC.       NOx.GT.        PT08.S3.NOx.       NO2.GT.     
 Min.   : 383.2   Min.   :   2.0   Min.   : 322.0   Min.   :  2.0  
 1st Qu.: 742.5   1st Qu.: 112.0   1st Qu.: 665.5   1st Qu.: 85.9  
 Median : 923.2   Median : 229.0   Median : 817.5   Median :113.1  
 Mean   : 939.0   Mean   : 246.9   Mean   : 835.4   Mean   :113.1  
 3rd Qu.:1104.8   3rd Qu.: 284.2   3rd Qu.: 960.2   3rd Qu.:133.0  
 Max.   :2214.0   Max.   :1479.0   Max.   :2682.8   Max.   :339.7  
  PT08.S4.NO2.   PT08.S5.O3.           T               RH        
 Min.   : 551   Min.   : 221.0   Min.   :-1.90   Min.   : 9.175  
 1st Qu.:1242   1st Qu.: 741.8   1st Qu.:12.03   1st Qu.:36.550  
 Median :1456   Median : 982.5   Median :18.27   Median :49.232  
 Mean   :1456   Mean   :1022.8   Mean   :18.32   Mean   :49.232  
 3rd Qu.:1662   3rd Qu.:1255.2   3rd Qu.:24.07   3rd Qu.:61.875  
 Max.   :2775   Max.   :2522.8   Max.   :44.60   Max.   :88.725  
       AH        
 Min.   :0.1847  
 1st Qu.:0.7461  
 Median :1.0154  
 Mean   :1.0255  
 3rd Qu.:1.2962  
 Max.   :2.2310

# Converting data to long format and scaling the values
air_quality_long <- air_quality %>%
  pivot_longer(cols = -c(Date, Time), names_to = "variable", values_to = "value") %>%
  mutate(value = scale(value))  # Scale the data

# 1. Histograms of Scaled Data
ggplot(air_quality_long, aes(x = value)) +
  geom_histogram(bins = 30, fill = "blue", color = "black") +
  facet_wrap(~ variable, scales = "free_x") +
  theme_minimal() +
  theme(plot.margin = margin(20, 20, 20, 20)) +
  labs(title = "Histograms of Scaled Air Quality Variables",
       x = "Scaled Value",
       y = "Frequency")

# 2. Box Plots of Scaled Data
ggplot(air_quality_long, aes(x = variable, y = value)) +
  geom_boxplot(fill = "blue", color = "black") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(title = "Box Plots of Scaled Air Quality Variables",
       x = "Variable",
       y = "Scaled Value")

names(air_quality)

ggplot(air_quality, aes(x = `CO.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of CO(GT)") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

Warning message:
“Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.”

ggplot(air_quality, aes(x = `CO.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of CO(GT)") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`CO.GT.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original CO.GT.", xlab = "CO.GT.", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed CO.GT.", xlab = "Transformed CO.GT.", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.343434343434343"

ggplot(air_quality, aes(x = `PT08.S1.CO.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of PT08.S1.CO.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `PT08.S1.CO.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale PT08.S1(CO)") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`PT08.S1.CO.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `PT08.S1.CO.`", xlab = "`PT08.S1.CO.`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `PT08.S1.CO.`", xlab = "Transformed `PT08.S1.CO.`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  -0.626262626262626"

ggplot(air_quality, aes(x = `NMHC.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of NMHC.GT.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `NMHC.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale NMHC.GT.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`NMHC.GT.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `NMHC.GT.`", xlab = "`NMHC.GT.`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `NMHC.GT.`", xlab = "Transformed `NMHC.GT.`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.585858585858586"

ggplot(air_quality, aes(x = `C6H6.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of C6H6.GT.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `C6H6.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale C6H6.GT.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`C6H6.GT.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `C6H6.GT.`", xlab = "`C6H6.GT.`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `C6H6.GT.`", xlab = "Transformed `C6H6.GT.`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.303030303030303"

ggplot(air_quality, aes(x = `PT08.S2.NMHC.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of PT08.S2.NMHC.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `PT08.S2.NMHC.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale PT08.S2.NMHC.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`PT08.S2.NMHC.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `PT08.S2.NMHC.`", xlab = "`PT08.S2.NMHC.`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `PT08.S2.NMHC.`", xlab = "Transformed `PT08.S2.NMHC.`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.181818181818182"

ggplot(air_quality, aes(x = `NOx.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of NOx.GT.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `NOx.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale NOx.GT.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`NOx.GT.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `NOx.GT.`", xlab = "`NOx.GT.`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `NOx.GT.`", xlab = "Transformed `NOx.GT.`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.222222222222222"

ggplot(air_quality, aes(x = `PT08.S3.NOx.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of PT08.S3.NOx.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `PT08.S3.NOx.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale PT08.S3.NOx.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`PT08.S3.NOx.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `PT08.S3.NOx.`", xlab = "`PT08.S3.NOx.`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `PT08.S3.NOx.`", xlab = "Transformed `PT08.S3.NOx.`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  -0.0202020202020201"

ggplot(air_quality, aes(x = `NO2.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of NO2.GT.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `NO2.GT.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale NO2.GT.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`NO2.GT.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `NO2.GT.`", xlab = "`NO2.GT.`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `NO2.GT.`", xlab = "Transformed `NO2.GT.`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.585858585858586"

ggplot(air_quality, aes(x = `PT08.S4.NO2.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of PT08.S4.NO2.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `PT08.S4.NO2.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale PT08.S4.NO2.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`PT08.S4.NO2.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `PT08.S4.NO2.`", xlab = "`PT08.S4.NO2.`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `PT08.S4.NO2.`", xlab = "Transformed `PT08.S4.NO2.`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.707070707070707"

ggplot(air_quality, aes(x = `PT08.S5.O3.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of PT08.S5.O3.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `PT08.S5.O3.`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale PT08.S5.O3.") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`PT08.S5.O3.`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `PT08.S5.O3.`", xlab = "`PT08.S5.O3.`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `PT08.S5.O3.`", xlab = "Transformed `PT08.S5.O3.`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.343434343434343"

ggplot(air_quality, aes(x = `T`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of T") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `T`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale T") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

Warning message in transformation$transform(x):
“NaNs produced”
Warning message in scale_x_log10():
“log-10 transformation introduced infinite values.”
Warning message in transformation$transform(x):
“NaNs produced”
Warning message in scale_x_log10():
“log-10 transformation introduced infinite values.”
Warning message:
“Removed 14 rows containing non-finite outside the scale range (`stat_bin()`).”
Warning message:
“Removed 14 rows containing non-finite outside the scale range
(`stat_density()`).”

response_variable <- air_quality$`T`

shift_value <- abs(min(response_variable)) + 1
response_variable <- response_variable + shift_value

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `T`", xlab = "`T`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `T`", xlab = "Transformed `T`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.626262626262626"

ggplot(air_quality, aes(x = `RH`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  theme_minimal() +
  ggtitle("Distribution of RH") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `RH`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale RH") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`RH`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `RH`", xlab = "`RH`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `RH`", xlab = "Transformed `RH`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.909090909090909"

ggplot(air_quality, aes(x = `AH`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "red", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of AH") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

ggplot(air_quality, aes(x = `AH`)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "cyan", color = "black", alpha = 0.6) +
  geom_density(color = "blue", size = 1) +
  scale_x_log10() +  # Apply log scale to the x-axis
  theme_minimal() +
  ggtitle("Distribution of logscale AH") +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 12),  # Increase font size for x-axis labels
    axis.text.y = element_text(size = 12)   # Increase font size for y-axis labels
  )

response_variable <- air_quality$`AH`

response_variable <- ifelse(response_variable == 0, response_variable + 1, response_variable)

boxcox_result <- boxcox(response_variable ~ 1, plotit = TRUE)

optimal_lambda <- boxcox_result$x[which.max(boxcox_result$y)]
print(paste("Optimal lambda: ", optimal_lambda))
if (optimal_lambda == 0) {
  transformed_response <- log(response_variable)
} else {
  transformed_response <- (response_variable^optimal_lambda - 1) / optimal_lambda
}

par(mfrow = c(1, 2)) # Split the plotting window into 1 row and 2 columns
hist(response_variable, main = "Original `AH`", xlab = "`AH`", ylab = "Frequency", col = "skyblue", border = "white")
hist(transformed_response, main = "Transformed `AH`", xlab = "Transformed `AH`", ylab = "Frequency", col = "skyblue", border = "white")
par(mfrow = c(1, 1))

[1] "Optimal lambda:  0.666666666666667"

cor_matrix <- cor(air_quality[, -c(14, 15)])
cor_matrix

library(corrplot)
corrplot(round(cor(cor_matrix), 2), method="number")

corrplot 0.92 loaded


Attaching package: ‘corrplot’


The following object is masked from ‘package:pls’:

    corrplot

threshold <- 0.7
high_cor_pairs <- list()

for(i in 1:(nrow(cor_matrix)-1)) {
  for(j in (i+1):ncol(cor_matrix)) {
    if(abs(cor_matrix[i, j]) > threshold) {
      high_cor_pairs <- rbind(high_cor_pairs, 
                              data.frame(Var1 = rownames(cor_matrix)[i], 
                                         Var2 = colnames(cor_matrix)[j], 
                                         value = cor_matrix[i, j]))
    }
  }
}

print(high_cor_pairs)

            Var1          Var2      value
1         CO.GT.   PT08.S1.CO.  0.7733943
2         CO.GT.      C6H6.GT.  0.8123923
3         CO.GT. PT08.S2.NMHC.  0.7955864
4         CO.GT.       NOx.GT.  0.7622972
5         CO.GT.   PT08.S5.O3.  0.7590266
6    PT08.S1.CO.      C6H6.GT.  0.8838209
7    PT08.S1.CO. PT08.S2.NMHC.  0.8929724
8    PT08.S1.CO.  PT08.S3.NOx. -0.7719181
9    PT08.S1.CO.   PT08.S5.O3.  0.8993259
10      C6H6.GT. PT08.S2.NMHC.  0.9819620
11      C6H6.GT.  PT08.S3.NOx. -0.7357112
12      C6H6.GT.  PT08.S4.NO2.  0.7657168
13      C6H6.GT.   PT08.S5.O3.  0.8657266
14 PT08.S2.NMHC.  PT08.S3.NOx. -0.7966873
15 PT08.S2.NMHC.  PT08.S4.NO2.  0.7772348
16 PT08.S2.NMHC.   PT08.S5.O3.  0.8805903
17       NOx.GT.       NO2.GT.  0.7631329
18  PT08.S3.NOx.   PT08.S5.O3. -0.7965536

threshold <- 0

neg_cor_pairs <- list()

for(i in 1:(nrow(cor_matrix)-1)) {
  for(j in (i+1):ncol(cor_matrix)) {
    if(cor_matrix[i, j] < threshold) {
      neg_cor_pairs <- rbind(neg_cor_pairs, 
                             data.frame(Var1 = rownames(cor_matrix)[i], 
                                        Var2 = colnames(cor_matrix)[j], 
                                        value = cor_matrix[i, j]))
    }
  }
}

print(neg_cor_pairs)

            Var1         Var2       value
1         CO.GT. PT08.S3.NOx. -0.61387028
2    PT08.S1.CO. PT08.S3.NOx. -0.77191812
3       NMHC.GT. PT08.S3.NOx. -0.26198519
4       NMHC.GT.           RH -0.05279411
5       C6H6.GT. PT08.S3.NOx. -0.73571121
6       C6H6.GT.           RH -0.06164347
7  PT08.S2.NMHC. PT08.S3.NOx. -0.79668731
8  PT08.S2.NMHC.           RH -0.09035172
9        NOx.GT. PT08.S3.NOx. -0.56325927
10       NOx.GT.            T -0.23565652
11       NOx.GT.           AH -0.12683093
12  PT08.S3.NOx.      NO2.GT. -0.56953458
13  PT08.S3.NOx. PT08.S4.NO2. -0.53846012
14  PT08.S3.NOx.  PT08.S5.O3. -0.79655364
15  PT08.S3.NOx.            T -0.14513301
16  PT08.S3.NOx.           RH -0.05672984
17  PT08.S3.NOx.           AH -0.23202063
18       NO2.GT.            T -0.16531710
19       NO2.GT.           RH -0.08064472
20       NO2.GT.           AH -0.29119971
21  PT08.S4.NO2.           RH -0.03218809
22   PT08.S5.O3.            T -0.02719336
23             T           RH -0.57856879

# Identify outliers using the IQR method
Q1 <- quantile(air_quality$'C6H6(GT)', 0.25)
Q3 <- quantile(air_quality$'C6H6(GT)', 0.75)
IQR <- Q3 - Q1
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
outliers <- which(air_quality$'C6H6(GT)' < lower_bound | air_quality$'C6H6(GT)' > upper_bound)

print(length(outliers))

lower_bound <- Q1 - 3 * IQR
upper_bound <- Q3 + 3 * IQR
extreme_outliers <- which(air_quality$'C6H6(GT)' < lower_bound | air_quality$'C6H6(GT)' > upper_bound)
print(length(extreme_outliers))

[1] 0
[1] 0

summary(air_quality)

     CO.GT.        PT08.S1.CO.        NMHC.GT.         C6H6.GT.     
 Min.   : 0.100   Min.   : 647.2   Min.   :   7.0   Min.   : 0.149  
 1st Qu.: 1.200   1st Qu.: 941.2   1st Qu.: 218.8   1st Qu.: 4.591  
 Median : 2.153   Median :1074.5   Median : 218.8   Median : 8.593  
 Mean   : 2.153   Mean   :1099.7   Mean   : 218.8   Mean   :10.083  
 3rd Qu.: 2.600   3rd Qu.:1221.2   3rd Qu.: 218.8   3rd Qu.:13.636  
 Max.   :11.900   Max.   :2039.8   Max.   :1189.0   Max.   :63.741  
 PT08.S2.NMHC.       NOx.GT.        PT08.S3.NOx.       NO2.GT.     
 Min.   : 383.2   Min.   :   2.0   Min.   : 322.0   Min.   :  2.0  
 1st Qu.: 742.5   1st Qu.: 112.0   1st Qu.: 665.5   1st Qu.: 85.9  
 Median : 923.2   Median : 229.0   Median : 817.5   Median :113.1  
 Mean   : 939.0   Mean   : 246.9   Mean   : 835.4   Mean   :113.1  
 3rd Qu.:1104.8   3rd Qu.: 284.2   3rd Qu.: 960.2   3rd Qu.:133.0  
 Max.   :2214.0   Max.   :1479.0   Max.   :2682.8   Max.   :339.7  
  PT08.S4.NO2.   PT08.S5.O3.           T               RH        
 Min.   : 551   Min.   : 221.0   Min.   :-1.90   Min.   : 9.175  
 1st Qu.:1242   1st Qu.: 741.8   1st Qu.:12.03   1st Qu.:36.550  
 Median :1456   Median : 982.5   Median :18.27   Median :49.232  
 Mean   :1456   Mean   :1022.8   Mean   :18.32   Mean   :49.232  
 3rd Qu.:1662   3rd Qu.:1255.2   3rd Qu.:24.07   3rd Qu.:61.875  
 Max.   :2775   Max.   :2522.8   Max.   :44.60   Max.   :88.725  
       AH              Date                       
 Min.   :0.1847   Min.   :2004-03-10 00:00:00.00  
 1st Qu.:0.7461   1st Qu.:2004-06-16 00:00:00.00  
 Median :1.0154   Median :2004-09-21 00:00:00.00  
 Mean   :1.0255   Mean   :2004-09-21 04:30:05.19  
 3rd Qu.:1.2962   3rd Qu.:2004-12-28 00:00:00.00  
 Max.   :2.2310   Max.   :2005-04-04 00:00:00.00  
      Time                       
 Min.   :1899-12-31 00:00:00.00  
 1st Qu.:1899-12-31 05:00:00.00  
 Median :1899-12-31 11:00:00.00  
 Mean   :1899-12-31 11:29:54.80  
 3rd Qu.:1899-12-31 18:00:00.00  
 Max.   :1899-12-31 23:00:00.00

tmp <- bc_transform(air_quality[, -c(14, 15)])
air_quality <- tmp

[1] "CO.GT."
[1] "Optimal lambda: 0.7"
[1] "PT08.S1.CO."
[1] "Optimal lambda: 0.7"
[1] "NMHC.GT."
[1] "Optimal lambda: 0.7"
[1] "C6H6.GT."
[1] "Optimal lambda: 0.7"
[1] "PT08.S2.NMHC."
[1] "Optimal lambda: 0.7"
[1] "NOx.GT."
[1] "Optimal lambda: 0.7"
[1] "PT08.S3.NOx."
[1] "Optimal lambda: 0.7"
[1] "NO2.GT."
[1] "Optimal lambda: 0.7"
[1] "PT08.S4.NO2."
[1] "Optimal lambda: 0.7"
[1] "PT08.S5.O3."
[1] "Optimal lambda: 0.7"
[1] "T"
[1] "Optimal lambda: 0.7"
[1] "RH"
[1] "Optimal lambda: 0.7"
[1] "AH"
[1] "Optimal lambda: 0.7"

split_ratio <- 0.7
split_index <- floor(nrow(air_quality) * split_ratio)

train = air_quality[1:split_index,]
test = air_quality[(split_index + 1):nrow(air_quality),]

# số chiều tập train
dim(train)
rownames(train) <- 1:nrow(train)

# số chiều tập test
dim(test)
rownames(test) <- 1:nrow(test)

# xem một số quan trắc của tập train
str(train)

'data.frame':	6549 obs. of  13 variables:
 $ CO.GT.       : num  1.36 0.892 1.052 1.052 0.557 ...
 $ PT08.S1.CO.  : num  222 214 226 223 211 ...
 $ NMHC.GT.     : num  46.2 37.4 31.4 29.3 21 ...
 $ C6H6.GT.     : num  6.65 5.43 5.22 5.34 3.88 ...
 $ PT08.S2.NMHC.: num  184 173 171 172 157 ...
 $ NOx.GT.      : num  49.7 35.2 41.9 51 41.9 ...
 $ PT08.S3.NOx. : num  185 200 196 190 203 ...
 $ NO2.GT.      : num  37.7 32.4 37.9 39.8 38.4 ...
 $ PT08.S4.NO2. : num  258 244 243 247 236 ...
 $ PT08.S5.O3.  : num  211 175 188 203 192 ...
 $ T            : num  8.74 8.61 7.99 7.59 7.66 ...
 $ RH           : num  20.3 19.9 21.9 23.7 23.5 ...
 $ AH           : num  -0.252 -0.287 -0.26 -0.221 -0.219 ...

full.lm <- lm(`C6H6.GT.` ~ ., data = train)
print(summary(full.lm))

Call:
lm(formula = C6H6.GT. ~ ., data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3855 -0.1963 -0.0553  0.1558  3.9380 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.423e+01  1.795e-01 -79.253  < 2e-16 ***
CO.GT.         4.689e-02  8.772e-03   5.346 9.32e-08 ***
PT08.S1.CO.   -1.518e-03  4.185e-04  -3.626 0.000290 ***
NMHC.GT.       1.064e-03  2.993e-04   3.556 0.000379 ***
PT08.S2.NMHC.  1.083e-01  5.774e-04 187.655  < 2e-16 ***
NOx.GT.        8.187e-03  2.838e-04  28.845  < 2e-16 ***
PT08.S3.NOx.   1.125e-02  3.328e-04  33.819  < 2e-16 ***
NO2.GT.       -1.848e-02  7.815e-04 -23.649  < 2e-16 ***
PT08.S4.NO2.   5.558e-03  3.845e-04  14.453  < 2e-16 ***
PT08.S5.O3.    1.357e-03  2.494e-04   5.443 5.44e-08 ***
T             -8.847e-02  7.468e-03 -11.846  < 2e-16 ***
RH            -4.064e-02  3.545e-03 -11.464  < 2e-16 ***
AH             5.389e-01  4.896e-02  11.008  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3119 on 6536 degrees of freedom
Multiple R-squared:  0.9923,	Adjusted R-squared:  0.9923 
F-statistic: 7.05e+04 on 12 and 6536 DF,  p-value: < 2.2e-16

# Mô hình chặn dưới
model.lb <- lm(`C6H6.GT.` ~ 1, data = train)

# Mô hình chặn trên
model.up <- full.lm

step(full.lm, scope = list(lower = model.lb, upper = model.up), direction = "both", trace = FALSE)

Call:
lm(formula = C6H6.GT. ~ CO.GT. + PT08.S1.CO. + NMHC.GT. + PT08.S2.NMHC. + 
    NOx.GT. + PT08.S3.NOx. + NO2.GT. + PT08.S4.NO2. + PT08.S5.O3. + 
    T + RH + AH, data = train)

Coefficients:
  (Intercept)         CO.GT.    PT08.S1.CO.       NMHC.GT.  PT08.S2.NMHC.  
   -14.229708       0.046891      -0.001518       0.001064       0.108347  
      NOx.GT.   PT08.S3.NOx.        NO2.GT.   PT08.S4.NO2.    PT08.S5.O3.  
     0.008187       0.011255      -0.018483       0.005558       0.001357  
            T             RH             AH  
    -0.088472      -0.040642       0.538945

air_models<- regsubsets(C6H6.GT. ~ CO.GT. + PT08.S1.CO. + NMHC.GT. + PT08.S2.NMHC. + 
    NOx.GT. + PT08.S3.NOx. + NO2.GT. + PT08.S4.NO2. + PT08.S5.O3. + 
    T + RH + AH, data = train)
summary.air<-summary(air_models)

# Lựa chọn mô hình tốt nhất từ reg subsets 
summary.air$which

summary.air$rsq

# model with largest adjusted R^2 
summary.air$adjr2

# model with smallest Mallow's Cp
summary.air$cp

# Tiêu chí chọn mô hình tốt nhất 4: mô hình với BIC nhỏ
summary.air$bic

best_model_index <- which.min(summary.air$bic)
best_model <- summary.air$which[best_model_index, ]
print(best_model)
best_vars <- names(best_model[best_model])
best_vars <- best_vars[best_vars != "(Intercept)"]
print(best_vars)

  (Intercept)        CO.GT.   PT08.S1.CO.      NMHC.GT. PT08.S2.NMHC. 
         TRUE         FALSE         FALSE         FALSE          TRUE 
      NOx.GT.  PT08.S3.NOx.       NO2.GT.  PT08.S4.NO2.   PT08.S5.O3. 
         TRUE          TRUE          TRUE          TRUE         FALSE 
            T            RH            AH 
         TRUE          TRUE          TRUE 
[1] "PT08.S2.NMHC." "NOx.GT."       "PT08.S3.NOx."  "NO2.GT."      
[5] "PT08.S4.NO2."  "T"             "RH"            "AH"

# Xây dựng mô hình tốt nhất
formula_str <- paste("C6H6.GT. ~", paste(best_vars, collapse = " + "))
best_model_air <- lm(as.formula(formula_str), data=train)

# Tóm tắt mô hình
summary(best_model_air)

Call:
lm(formula = as.formula(formula_str), data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3449 -0.1997 -0.0546  0.1609  3.8753 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.434e+01  1.717e-01  -83.53   <2e-16 ***
PT08.S2.NMHC.  1.091e-01  5.156e-04  211.67   <2e-16 ***
NOx.GT.        8.707e-03  2.648e-04   32.88   <2e-16 ***
PT08.S3.NOx.   1.081e-02  3.248e-04   33.27   <2e-16 ***
NO2.GT.       -1.743e-02  7.717e-04  -22.59   <2e-16 ***
PT08.S4.NO2.   5.917e-03  3.628e-04   16.31   <2e-16 ***
T             -9.191e-02  7.381e-03  -12.45   <2e-16 ***
RH            -4.096e-02  3.558e-03  -11.51   <2e-16 ***
AH             5.359e-01  4.887e-02   10.97   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3136 on 6540 degrees of freedom
Multiple R-squared:  0.9922,	Adjusted R-squared:  0.9922 
F-statistic: 1.046e+05 on 8 and 6540 DF,  p-value: < 2.2e-16

# Trực quan hóa
par(mfrow = c(2, 2))
plot(best_model_air)

plot(best_model_air, which=1, col=c("blue")) # Residuals vs Fitted Plot

plot(best_model_air, which=2, col=c("blue")) # QQ Plot

# Shapiro-Wilk normality test
# shapiro.test(residuals(best_model_air)[3:5000])
CheckNormal(best_model_air)

	Shapiro-Wilk normality test

data:  model$residuals[3:5000]
W = 0.96092, p-value < 2.2e-16

[1] "H0 rejected: the residuals are NOT distributed normally"

plot(best_model_air, which=3, col=c("blue")) # Scale-Location

# Breusch-Pagan Test
# bptest(best_model_air)
CheckHomos(best_model_air)

	studentized Breusch-Pagan test

data:  model
BP = 1157, df = 8, p-value < 2.2e-16

[1] "H0 rejected: Error variance spreads INCONSTANTLY/generating patterns (Heteroscedasticity)"

plot(best_model_air, which=5, col=c("blue")) # Residuals vs Leverage

plot(best_model_air, which=4, col=c("blue"))

# Xây dựng ngưỡng cho Cook Distance
threshold <- 4 / nrow(train)
threshold

# Tính toán Cook Distance
cooks_d <- cooks.distance(best_model_air)

# Xác định các influential point dựa trên ngưỡng
influential_points <- which(cooks_d > threshold)

# Xây dựng model
best_model_air.2 <- lm(as.formula(formula_str), data=train[-c(influential_points), ])
summary(best_model_air.2)

# Shapiro-Wilk normality test
CheckNormal(best_model_air.2)

# Breusch-Pagan Test
CheckHomos(best_model_air.2)

Call:
lm(formula = as.formula(formula_str), data = train[-c(influential_points), 
    ])

Residuals:
     Min       1Q   Median       3Q      Max 
-0.56389 -0.15826 -0.04263  0.13520  0.80613 

Coefficients:
                Estimate Std. Error  t value Pr(>|t|)    
(Intercept)   -1.481e+01  1.445e-01 -102.490  < 2e-16 ***
PT08.S2.NMHC.  1.086e-01  4.857e-04  223.638  < 2e-16 ***
NOx.GT.        6.541e-03  2.133e-04   30.671  < 2e-16 ***
PT08.S3.NOx.   9.981e-03  3.296e-04   30.287  < 2e-16 ***
NO2.GT.       -1.143e-02  5.982e-04  -19.113  < 2e-16 ***
PT08.S4.NO2.   3.640e-03  3.127e-04   11.639  < 2e-16 ***
T             -3.867e-02  5.673e-03   -6.816 1.02e-11 ***
RH            -1.452e-02  2.750e-03   -5.278 1.35e-07 ***
AH             3.109e-01  3.802e-02    8.178 3.49e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2266 on 6160 degrees of freedom
Multiple R-squared:  0.9947,	Adjusted R-squared:  0.9947 
F-statistic: 1.449e+05 on 8 and 6160 DF,  p-value: < 2.2e-16

	Shapiro-Wilk normality test

data:  model$residuals[3:5000]
W = 0.96912, p-value < 2.2e-16

[1] "H0 rejected: the residuals are NOT distributed normally"

	studentized Breusch-Pagan test

data:  model
BP = 801.16, df = 8, p-value < 2.2e-16

[1] "H0 rejected: Error variance spreads INCONSTANTLY/generating patterns (Heteroscedasticity)"

# Tính toán Cook Distance
cooks_d <- cooks.distance(best_model_air.2)

# Xác định các influential point dựa trên ngưỡng
influential_points <- which(cooks_d > threshold)

# Xây dựng model
best_model_air.3 <- lm(as.formula(formula_str), data=train[-c(influential_points), ])
summary(best_model_air.3)

# Shapiro-Wilk normality test
CheckNormal(best_model_air.3)

# Breusch-Pagan Test
CheckHomos(best_model_air.3)

Call:
lm(formula = as.formula(formula_str), data = train[-c(influential_points), 
    ])

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3360 -0.1997 -0.0562  0.1612  3.8857 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.434e+01  1.766e-01  -81.18   <2e-16 ***
PT08.S2.NMHC.  1.089e-01  5.315e-04  204.97   <2e-16 ***
NOx.GT.        8.682e-03  2.730e-04   31.80   <2e-16 ***
PT08.S3.NOx.   1.064e-02  3.348e-04   31.79   <2e-16 ***
NO2.GT.       -1.733e-02  8.020e-04  -21.61   <2e-16 ***
PT08.S4.NO2.   5.936e-03  3.744e-04   15.86   <2e-16 ***
T             -8.964e-02  7.607e-03  -11.78   <2e-16 ***
RH            -3.969e-02  3.665e-03  -10.83   <2e-16 ***
AH             5.084e-01  5.031e-02   10.10   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3133 on 6095 degrees of freedom
Multiple R-squared:  0.9922,	Adjusted R-squared:  0.9922 
F-statistic: 9.736e+04 on 8 and 6095 DF,  p-value: < 2.2e-16

	Shapiro-Wilk normality test

data:  model$residuals[3:5000]
W = 0.94602, p-value < 2.2e-16

[1] "H0 rejected: the residuals are NOT distributed normally"

	studentized Breusch-Pagan test

data:  model
BP = 1088.7, df = 8, p-value < 2.2e-16

[1] "H0 rejected: Error variance spreads INCONSTANTLY/generating patterns (Heteroscedasticity)"

# Tính toán Cook Distance
cooks_d <- cooks.distance(best_model_air.3)

# Xác định các influential point dựa trên ngưỡng
influential_points <- which(cooks_d > threshold)

# Xây dựng model
best_model_air.4 <- lm(as.formula(formula_str), data=train[-c(influential_points), ])
summary(best_model_air.4)

# Shapiro-Wilk normality test
CheckNormal(best_model_air.4)

# Breusch-Pagan Test
CheckHomos(best_model_air.4)

Call:
lm(formula = as.formula(formula_str), data = train[-c(influential_points), 
    ])

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3138 -0.2009 -0.0543  0.1611  3.9055 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.430e+01  1.748e-01  -81.82   <2e-16 ***
PT08.S2.NMHC.  1.087e-01  5.285e-04  205.62   <2e-16 ***
NOx.GT.        8.774e-03  2.735e-04   32.08   <2e-16 ***
PT08.S3.NOx.   1.044e-02  3.328e-04   31.36   <2e-16 ***
NO2.GT.       -1.730e-02  7.979e-04  -21.68   <2e-16 ***
PT08.S4.NO2.   6.157e-03  3.736e-04   16.48   <2e-16 ***
T             -8.998e-02  7.540e-03  -11.94   <2e-16 ***
RH            -4.046e-02  3.646e-03  -11.10   <2e-16 ***
AH             5.151e-01  5.005e-02   10.29   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3137 on 6157 degrees of freedom
Multiple R-squared:  0.9922,	Adjusted R-squared:  0.9922 
F-statistic: 9.837e+04 on 8 and 6157 DF,  p-value: < 2.2e-16

	Shapiro-Wilk normality test

data:  model$residuals[3:5000]
W = 0.95024, p-value < 2.2e-16

[1] "H0 rejected: the residuals are NOT distributed normally"

	studentized Breusch-Pagan test

data:  model
BP = 1042, df = 8, p-value < 2.2e-16

[1] "H0 rejected: Error variance spreads INCONSTANTLY/generating patterns (Heteroscedasticity)"

# Tính toán Cook Distance
cooks_d <- cooks.distance(best_model_air.4)

# Xác định các influential point dựa trên ngưỡng
influential_points <- which(cooks_d > threshold)

# Xây dựng model
best_model_air.5 <- lm(as.formula(formula_str), data=train[-c(influential_points), ])
summary(best_model_air.4)

# Shapiro-Wilk normality test
CheckNormal(best_model_air.4)

# Breusch-Pagan Test
CheckHomos(best_model_air.4)

Call:
lm(formula = as.formula(formula_str), data = train[-c(influential_points), 
    ])

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3138 -0.2009 -0.0543  0.1611  3.9055 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.430e+01  1.748e-01  -81.82   <2e-16 ***
PT08.S2.NMHC.  1.087e-01  5.285e-04  205.62   <2e-16 ***
NOx.GT.        8.774e-03  2.735e-04   32.08   <2e-16 ***
PT08.S3.NOx.   1.044e-02  3.328e-04   31.36   <2e-16 ***
NO2.GT.       -1.730e-02  7.979e-04  -21.68   <2e-16 ***
PT08.S4.NO2.   6.157e-03  3.736e-04   16.48   <2e-16 ***
T             -8.998e-02  7.540e-03  -11.94   <2e-16 ***
RH            -4.046e-02  3.646e-03  -11.10   <2e-16 ***
AH             5.151e-01  5.005e-02   10.29   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3137 on 6157 degrees of freedom
Multiple R-squared:  0.9922,	Adjusted R-squared:  0.9922 
F-statistic: 9.837e+04 on 8 and 6157 DF,  p-value: < 2.2e-16

	Shapiro-Wilk normality test

data:  model$residuals[3:5000]
W = 0.95024, p-value < 2.2e-16

[1] "H0 rejected: the residuals are NOT distributed normally"

	studentized Breusch-Pagan test

data:  model
BP = 1042, df = 8, p-value < 2.2e-16

[1] "H0 rejected: Error variance spreads INCONSTANTLY/generating patterns (Heteroscedasticity)"

# Predict using the model and evaluate on the test set
predictions <- predict(best_model_air.5, test)

# Compare predictions with actual values
plot(test$`C6H6.GT.`, predictions, xlab = "Actual", ylab = "Predicted", main = "Predicted vs Actual C6H6(GT) Values (Multiple Regression Model)")  # Plot actual vs predicted values
abline(0, 1)  # Add a diagonal line for reference

# Calculate and print the Root Mean Squared Error (RMSE)
rmse <- sqrt(mean((test$`C6H6.GT.` - predictions)^2))  # Calculate RMSE between actual and predicted values
print(paste("RMSE: ", rmse))

[1] "RMSE:  0.520176239781131"

# Calculate the sum of squares of residuals
ss_res <- sum((test$`C6H6.GT.` - predictions)^2)

# Calculate the total sum of squares
ss_tot <- sum((test$`C6H6.GT.` - mean(test$`C6H6.GT.`))^2)

# Calculate R-squared
r_squared <- 1 - (ss_res / ss_tot)

# Print R-squared
print(paste("R-squared: ", r_squared))

[1] "R-squared:  0.974160108813221"

# Fitting the PCR model on the training data
pcr_model <- pcr(`C6H6.GT.` ~ ., data = train, scale = TRUE, validation = "CV") # Fit PCR model with cross-validation

summary(pcr_model)

Data: 	X dimension: 6549 12 
	Y dimension: 6549 1
Fit method: svdpc
Number of components considered: 12

VALIDATION: RMSEP
Cross-validated using 10 random segments.
       (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
CV           3.559    1.063    1.053    1.016    1.001   0.7721   0.7637
adjCV        3.559    1.063    1.053    1.016    1.001   0.7721   0.7637
       7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
CV      0.7315   0.7142   0.5740    0.5741    0.3266    0.3129
adjCV   0.7314   0.7141   0.5738    0.5740    0.3265    0.3128

TRAINING: % variance explained
          1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
X           51.80    68.66    80.52    87.53    92.94    95.34    97.08
C6H6.GT.    91.09    91.26    91.86    92.10    95.30    95.41    95.79
          8 comps  9 comps  10 comps  11 comps  12 comps
X           98.14    98.95     99.65     99.89    100.00
C6H6.GT.    95.99    97.41     97.41     99.17     99.23

# Plotting the RMSEP (Root Mean Squared Error of Prediction) to find the optimal number of components
validationplot(pcr_model, val.type = "RMSEP", main = "RMSEP (PCR Model)")

# Predict using the model and evaluate on the test set with optimal number of components
optimal_number_of_components <- 12  # Optimal number of components based on the RMSEP plot and summary
predictions <- predict(pcr_model, ncomp = optimal_number_of_components, newdata = test)

# Compare predictions with actual values
plot(test$`C6H6.GT.`, predictions, xlab = "Actual", ylab = "Predicted", main = "Predicted vs Actual C6H6(GT) Values (PCR Model)")  # Plot actual vs predicted values
abline(0, 1)  # Add a diagonal line for reference

# Calculate and print the Root Mean Squared Error (RMSE)
rmse <- sqrt(mean((test$`C6H6.GT.` - predictions)^2))  # Calculate RMSE between actual and predicted values
print(paste("RMSE: ", rmse))

[1] "RMSE:  0.517225894818149"

# Calculate the sum of squares of residuals
ss_res <- sum((test$`C6H6.GT.` - predictions)^2)

# Calculate the total sum of squares
ss_tot <- sum((test$`C6H6.GT.` - mean(test$`C6H6.GT.`))^2)

# Calculate R-squared
r_squared <- 1 - (ss_res / ss_tot)

# Print R-squared
print(paste("R-squared: ", r_squared))

[1] "R-squared:  0.97445239587654"

# Fitting the PLS model on the training data
pls_model <- plsr(`C6H6.GT.` ~ ., data = train, scale = TRUE, validation = "CV")

summary(pls_model)

Data: 	X dimension: 6549 12 
	Y dimension: 6549 1
Fit method: kernelpls
Number of components considered: 12

VALIDATION: RMSEP
Cross-validated using 10 random segments.
       (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
CV           3.559   0.9928   0.7452   0.6638   0.6035   0.4785   0.4408
adjCV        3.559   0.9928   0.7450   0.6638   0.6035   0.4783   0.4407
       7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
CV      0.3843   0.3322   0.3156    0.3145     0.313    0.3131
adjCV   0.3842   0.3321   0.3156    0.3144     0.313    0.3130

TRAINING: % variance explained
          1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
X           51.76    60.82    72.13    82.96    87.34    92.87    95.43
C6H6.GT.    92.23    95.63    96.54    97.14    98.21    98.48    98.84
          8 comps  9 comps  10 comps  11 comps  12 comps
X           96.64    98.12     99.15     99.30    100.00
C6H6.GT.    99.14    99.22     99.23     99.23     99.23

# Plotting the RMSEP (Root Mean Squared Error of Prediction) to find the optimal number of components
validationplot(pls_model, val.type = "RMSEP", main = "RMSEP (PLS Model)")

# Predict using the model and evaluate on the test set with optimal number of components
optimal_number_of_components <- 11  # Optimal number of components based on the RMSEP plot and summary
predictions2 <- predict(pls_model, ncomp = optimal_number_of_components, newdata = test)

# Compare predictions with actual values
plot(test$`C6H6.GT.`, predictions2, xlab = "Actual", ylab = "Predicted", main = "Predicted vs Actual C6H6(GT) Values (PLS Model)")  # Plot actual vs predicted values
abline(0, 1)  # Add a diagonal line for reference

# Calculate and print the Root Mean Squared Error (RMSE)
rmse2 <- sqrt(mean((test$`C6H6.GT.` - predictions2)^2))  # Calculate RMSE between actual and predicted values
print(paste("RMSE: ", rmse2))

[1] "RMSE:  0.517942615593671"

# Calculate the sum of squares of residuals
ss_res2 <- sum((test$`C6H6.GT.` - predictions2)^2)

# Calculate the total sum of squares
ss_tot2 <- sum((test$`C6H6.GT.` - mean(test$`C6H6.GT.`))^2)

# Calculate R-squared
r_squared2 <- 1 - (ss_res2 / ss_tot2)

# Print R-squared
print(paste("R-squared: ", r_squared2))

[1] "R-squared:  0.97438154410648"

library(randomForest)

set.seed(123)
model_rf <- randomForest(x = train[,-c(4)],
                         y = train$`C6H6.GT.`, 
                         ntree = 500)

model_rf

randomForest 4.7-1.1

Type rfNews() to see new features/changes/bug fixes.


Attaching package: ‘randomForest’


The following object is masked from ‘package:ggplot2’:

    margin


The following object is masked from ‘package:dplyr’:

    combine

Call:
 randomForest(x = train[, -c(4)], y = train$C6H6.GT., ntree = 500) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 4

          Mean of squared residuals: 0.02391061
                    % Var explained: 99.81

library(MLmetrics)
library(performance)

pred_rf_val <- predict(object = model_rf, newdata = test)

# Compare predictions with actual values
plot(test$`C6H6.GT.`, pred_rf_val, xlab = "Actual", ylab = "Predicted", main = "Predicted vs Actual C6H6(GT) Values (Random Forest)")  # Plot actual vs predicted values
abline(0, 1)  # Add a diagonal line for reference

# Calculate and print the Root Mean Squared Error (RMSE)
rmse2 <- sqrt(mean((test$`C6H6.GT.` - pred_rf_val)^2))  # Calculate RMSE between actual and predicted values
print(paste("RMSE: ", rmse2))

[1] "RMSE:  0.414118902326439"

# Calculate the sum of squares of residuals
ss_res2 <- sum((test$`C6H6.GT.` - pred_rf_val)^2)

# Calculate the total sum of squares
ss_tot2 <- sum((test$`C6H6.GT.` - mean(test$`C6H6.GT.`))^2)

# Calculate R-squared
r_squared2 <- 1 - (ss_res2 / ss_tot2)

# Print R-squared
print(paste("R-squared: ", r_squared2))

[1] "R-squared:  0.98362279569018"

library("tibble")
model_rf$importance %>% 
  as.data.frame() %>% 
  arrange(-IncNodePurity) %>% 
  rownames_to_column("variable") %>% 
  head(10) %>% 
  ggplot(aes(IncNodePurity, 
             reorder(variable, IncNodePurity))
         ) +
  geom_col(fill = "firebrick") +
  labs(x = "Importance",
       y = NULL,
       title = "Random Forest Variable Importance")

library(lime)

set.seed(123)
explainer <- lime(x = test[,-c(4)],
                  model = model_rf)

model_type.randomForest <- function(x){
  return("regression") # for regression problem
}

predict_model.randomForest <- function(x, newdata, type = "response") {

    # return prediction value
    predict(x, newdata) %>% as.data.frame()
    
}

#  Select only the first 4 observations
selected_data <- test[,-c(4)] %>% 
  slice(1:4)

#  Explain the model
set.seed(123)
explanation <- explain(x = selected_data, 
                       explainer = explainer,
                       n_features = 27 #  Number of features to explain the model
                       )

Warning message:
“Data contains numeric columns with zero variance”

Warning message in gower_work(x = x, y = y, pair_x = pair_x, pair_y = pair_y, n = NULL, :
“skipping variable with zero or non-finite range.”
Warning message in gower_work(x = x, y = y, pair_x = pair_x, pair_y = pair_y, n = NULL, :
“skipping variable with zero or non-finite range.”
Warning message in gower_work(x = x, y = y, pair_x = pair_x, pair_y = pair_y, n = NULL, :
“skipping variable with zero or non-finite range.”
Warning message in gower_work(x = x, y = y, pair_x = pair_x, pair_y = pair_y, n = NULL, :
“skipping variable with zero or non-finite range.”

library(plotly)

plot_features(explanation = explanation)

library(e1071)
model_svm <- svm(`C6H6.GT.` ~ ., data = train)
pred_svm_val <- predict(object = model_svm, newdata = test)

# Compare predictions with actual values
plot(test$`C6H6.GT.`, pred_svm_val, xlab = "Actual", ylab = "Predicted", main = "Predicted vs Actual C6H6(GT) Values (SVM)")  # Plot actual vs predicted values
abline(0, 1)  # Add a diagonal line for reference

# Calculate and print the Root Mean Squared Error (RMSE)
rmse2 <- sqrt(mean((test$`C6H6.GT.` - pred_svm_val)^2))  # Calculate RMSE between actual and predicted values
print(paste("RMSE: ", rmse2))

[1] "RMSE:  0.731968766359617"

# Calculate the sum of squares of residuals
ss_res2 <- sum((test$`C6H6.GT.` - pred_svm_val)^2)

# Calculate the total sum of squares
ss_tot2 <- sum((test$`C6H6.GT.` - mean(test$`C6H6.GT.`))^2)

# Calculate R-squared
r_squared2 <- 1 - (ss_res2 / ss_tot2)

# Print R-squared
print(paste("R-squared: ", r_squared2))

[1] "R-squared:  0.948834790325647"

library(lime)

# create the explanation for the SVR model.
set.seed(123)
explainer_svm <- lime(x = train[,-c(4)], 
                  model = model_svm)

# Create SVR model specification for lime.
model_type.svm <- function(x){
  return("regression") # for regression problem
}

predict_model.svm <- function(x, newdata, type = "response") {

    # return prediction value
    predict(x, newdata) %>% as.data.frame()
    
}

set.seed(123)
explanation_svm <- explain(x = selected_data, 
                       explainer = explainer_svm,
                       kernel_width = 1,
                       feature_select = "auto", # Method of feature selection for lime
                       n_features = 10 # Number of features to explain the model
                       )

set.seed(123)
explanation_svm <- explain(x = selected_data, 
                       explainer = explainer_svm,
                       kernel_width = 1,
                       feature_select = "auto", # Method of feature selection for lime
                       n_features = 10 # Number of features to explain the model
                       )

plot_features(explanation_svm)

Tên biến	Vai trò	Phân loại	Mô tả	Đơn vị	Missing Values
Date	Feature	Date (Ngày tháng năm)	Ngày tháng mà các giá trị độ đo được thu thập.		Không
Time	Feature	Categorical (Phân loại)	Thời gian mà các giá trị độ đo được thu thập.		Không
CO(GT)	Feature	Integer (Nguyên)	Nồng độ CO (carbon monoxide) trung bình thực tế theo giờ tính bằng mg/m^3 (máy phân tích tham chiếu)	mg/m^3	Không
PT08.S1(CO)	Feature	Categorical (Phân loại)	(Trung bình theo giờ) Phản hồi cảm biến (PPM) của cảm biến hồng ngoại không phân tán (NDIR) đối với carbon monoxide.		Không
NMHC(GT)	Feature	Integer (Nguyên)	Nồng độ Non Metanic HydroCarbons trong không khí (đơn vị microg/m^3/ trung bình theo giờ/ dựa trên máy phân tích tham chiếu)	microg/m^3	Không
C6H6(GT)	Feature	Continuous	Nồng độ trung bình Benzene trong không khí (đơn vị microg/m^3/ trung bình theo giờ/ dựa trên máy phân tích tham chiếu)	microg/m^3	Không
PT08.S2(NMHC)	Feature	Categorical (Phân loại)	(Trung bình theo giờ) Phản hồi cảm biến (PPM) của cảm biến hồng ngoại không phân tán (NDIR) đối với Non Metanic HydroCarbons		Không
NOx(GT)	Feature	Integer (Nguyên)	Nồng độ trung bình NOx (oxit nitơ) trong không khí (đơn vị PPB/ trung bình theo giờ/ dựa trên máy phân tích tham chiếu)	ppb	Không
PT08.S3(NOx)	Feature	Categorical (Phân loại)	(Trung bình theo giờ) Phản hồi cảm biến (PPM) của cảm biến hồng ngoại không phân tán (NDIR) đối với NOx (oxit nitơ)		Không
NO2(GT)	Feature	Integer (Nguyên)	Nồng độ trung bình NO2 (nitơ dioxit) trong không khí (đơn vị microg/m^3 / trung bình theo giờ/ dựa trên máy phân tích tham chiếu)	microg/m^3	Không
PT08.S4(NO2)	Feature	Categorical (Phân loại)	(Trung bình theo giờ) Phản hồi cảm biến (PPM) của cảm biến hồng ngoại không phân tán (NDIR) đối với NO2 (nitơ dioxit)		Không
PT08.S5(O3)	Feature	Categorical (Phân loại)	(Trung bình theo giờ) Phản hồi cảm biến (PPM) của cảm biến hồng ngoại không phân tán (NDIR) đối với O3		Không
T	Feature	Continuous (Liên tục)	Nhiệt độ (Thang độ Celsius)	°C	Không
RH	Feature	Continuous (Liên tục)	Độ ẩm tương đối (Relative Humidity)	%	Không
AH	Feature	Continuous (Liên tục)	Độ ẩm tuyệt đối (Absolute Humidity)		Không

	CO.GT.	PT08.S1.CO.	NMHC.GT.	C6H6.GT.	PT08.S2.NMHC.	NOx.GT.	PT08.S3.NOx.	NO2.GT.	PT08.S4.NO2.	PT08.S5.O3.	T	RH	AH
CO.GT.	1.00000000	0.7733943	0.28997516	0.81239231	0.79558637	0.76229720	-0.61387028	0.64652757	0.55822459	0.75902658	0.02026043	0.04320350	0.04327259
PT08.S1.CO.	0.77339429	1.0000000	0.27828805	0.88382089	0.89297245	0.61597363	-0.77191812	0.56299691	0.68287389	0.89932594	0.04862760	0.11460732	0.13532268
NMHC.GT.	0.28997516	0.2782880	1.00000000	0.28192285	0.27975111	0.09818457	-0.26198519	0.15654891	0.23298307	0.24059378	0.06680899	-0.05279411	0.03717300
C6H6.GT.	0.81239231	0.8838209	0.28192285	1.00000000	0.98196203	0.61687961	-0.73571121	0.53617833	0.76571684	0.86572663	0.19889091	-0.06164347	0.16790885
PT08.S2.NMHC.	0.79558637	0.8929724	0.27975111	0.98196203	1.00000000	0.60181206	-0.79668731	0.56142083	0.77723483	0.88059027	0.24132968	-0.09035172	0.18690569
NOx.GT.	0.76229720	0.6159736	0.09818457	0.61687961	0.60181206	1.00000000	-0.56325927	0.76313286	0.20552565	0.68821380	-0.23565652	0.19257781	-0.12683093
PT08.S3.NOx.	-0.61387028	-0.7719181	-0.26198519	-0.73571121	-0.79668731	-0.56325927	1.00000000	-0.56953458	-0.53846012	-0.79655364	-0.14513301	-0.05672984	-0.23202063
NO2.GT.	0.64652757	0.5629969	0.15654891	0.53617833	0.56142083	0.76313286	-0.56953458	1.00000000	0.14093950	0.62956418	-0.16531710	-0.08064472	-0.29119971
PT08.S4.NO2.	0.55822459	0.6828739	0.23298307	0.76571684	0.77723483	0.20552565	-0.53846012	0.14093950	1.00000000	0.59113739	0.56133337	-0.03218809	0.62964357
PT08.S5.O3.	0.75902658	0.8993259	0.24059378	0.86572663	0.88059027	0.68821380	-0.79655364	0.62956418	0.59113739	1.00000000	-0.02719336	0.12495163	0.07073801
T	0.02026043	0.0486276	0.06680899	0.19889091	0.24132968	-0.23565652	-0.14513301	-0.16531710	0.56133337	-0.02719336	1.00000000	-0.57856879	0.65647138
RH	0.04320350	0.1146073	-0.05279411	-0.06164347	-0.09035172	0.19257781	-0.05672984	-0.08064472	-0.03218809	0.12495163	-0.57856879	1.00000000	0.16797092
AH	0.04327259	0.1353227	0.03717300	0.16790885	0.18690569	-0.12683093	-0.23202063	-0.29119971	0.62964357	0.07073801	0.65647138	0.16797092	1.00000000

Date	Time	CO(GT)	PT08.S1(CO)	NMHC(GT)	C6H6(GT)	PT08.S2(NMHC)	NOx(GT)	PT08.S3(NOx)	NO2(GT)	PT08.S4(NO2)	PT08.S5(O3)	T	RH	AH
<dttm>	<dttm>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
2004-03-10	1899-12-31 18:00:00	2.6	1360.00	150	11.881723	1045.50	166	1056.25	113	1692.00	1267.50	13.600	48.875	0.7577538
2004-03-10	1899-12-31 19:00:00	2.0	1292.25	112	9.397165	954.75	103	1173.75	92	1558.75	972.25	13.300	47.700	0.7254874
2004-03-10	1899-12-31 20:00:00	2.2	1402.00	88	8.997817	939.25	131	1140.00	114	1554.50	1074.00	11.900	53.975	0.7502391
2004-03-10	1899-12-31 21:00:00	2.2	1375.50	80	9.228796	948.25	172	1092.00	122	1583.75	1203.25	11.000	60.000	0.7867125
2004-03-10	1899-12-31 22:00:00	1.6	1272.25	51	6.518224	835.50	131	1205.00	116	1490.00	1110.00	11.150	59.575	0.7887942
2004-03-10	1899-12-31 23:00:00	1.2	1197.00	38	4.741012	750.25	89	1336.50	96	1393.00	949.25	11.175	59.175	0.7847717

	(Intercept)	CO.GT.	PT08.S1.CO.	NMHC.GT.	PT08.S2.NMHC.	NOx.GT.	PT08.S3.NOx.	NO2.GT.	PT08.S4.NO2.	PT08.S5.O3.	T	RH	AH
1	TRUE	FALSE	FALSE	FALSE	TRUE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE
2	TRUE	FALSE	FALSE	FALSE	TRUE	FALSE	TRUE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE
3	TRUE	FALSE	FALSE	FALSE	TRUE	FALSE	TRUE	FALSE	FALSE	FALSE	FALSE	TRUE	FALSE
4	TRUE	FALSE	FALSE	FALSE	TRUE	TRUE	TRUE	TRUE	FALSE	FALSE	FALSE	FALSE	FALSE
5	TRUE	FALSE	FALSE	FALSE	TRUE	TRUE	TRUE	TRUE	TRUE	FALSE	FALSE	FALSE	FALSE
6	TRUE	FALSE	FALSE	FALSE	TRUE	TRUE	TRUE	TRUE	TRUE	TRUE	FALSE	FALSE	FALSE
7	TRUE	TRUE	FALSE	FALSE	TRUE	TRUE	TRUE	TRUE	TRUE	TRUE	FALSE	FALSE	FALSE
8	TRUE	FALSE	FALSE	FALSE	TRUE	TRUE	TRUE	TRUE	TRUE	FALSE	TRUE	TRUE	TRUE

PHÂN TÍCH CHẤT LƯỢNG KHÔNG KHÍ¶

Giới thiệu chung¶

Phát biểu bài toán¶

Giới thiệu về dữ liệu¶

Import thư viện¶

Hàm phụ trợ¶

Đọc dữ liệu¶

Khám phá và tiền xử lý dữ liệu¶

Dữ liệu có bao nhiều dòng và bao nhiêu cột?¶

Mỗi dòng có ý nghĩa gì? Liệu có tồn tại dòng nào mà mang ý nghĩa khác với các dòng còn lại không?¶

Dữ liệu có bị trùng lặp không?¶

Mỗi cột mang ý nghĩa gì?¶

Có cột nào có kiểu dữ liệu chưa phù hợp? Nếu có, cần chuyển đổi sang kiểu phù hợp¶

Các cột với kiểu dữ liệu số phân bố như thế nào?¶

Xử lý các giá trị bị thiếu¶

Phân tích đơn biến¶

1. CO(GT): Carbon monoxide concentration (mg/m³)¶

2. PT08.S1(CO): Sensor response for CO¶

3. NMHC(GT): Non-methane hydrocarbons concentration (μg/m³)¶

4. C6H6(GT): Benzene concentration (μg/m³)¶

5. PT08.S2(NMHC): Sensor response for NMHC¶

6. NOx(GT): Nitrogen oxides concentration (ppb)¶

7. PT08.S3(NOx): Sensor response for NOx¶

8. NO2(GT): Nitrogen dioxide concentration (μg/m³)¶

9. PT08.S4(NO2): Sensor response for NO2¶

10. PT08.S5(O3): Sensor response for ozone¶

11. T: Temperature (°C)¶

12. RH: Relative humidity (%)¶

13. AH: Absolute humidity (g/m³)¶

Phân tích đa biến¶

Tương quan giữa các biến¶

Khảo sát ngoại lai¶

Phân chia tập dữ liệu¶

Mô hình hóa¶

Mô hình đầy đủ¶

Lựa chọn mô hình¶

Phân tích Residuals vs Fitted Plot¶

Phân tích Normal Q–Q (quantile-quantile) Plot¶

Phân tích Scale-Location¶

Phân tích Residuals vs Leverage¶

Loại bỏ các điểm ảnh hưởng dựa trên Cook Distance¶

Dự đoán và đánh giá kết quả¶

Mô hình hóa bằng PCR¶

Mô hình hóa bằng PLS¶

Random forest¶

Support Vector Machine¶

So sánh các mô hình¶

Kết luận¶